% In this practical we will work with a single subject's EEG dat and perform an ERF
% analysis in sensor space. 
%

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

clear all; close all;

global OSLDIR;


tilde='/Users/dantemant/Documents/meeg';
osldir=[tilde '/osl1.2.beta.17'];    

tilde='/Users/woolrich/homedir';
osldir=[tilde '/matlab/osl1.2.beta.17'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

%testdir=[tilde '/meeg_data/meeg_signals/case_1084/111202'];
%testdir=['/Users/dantemant/Documents/meeg/duncan_data'];
testdir=[tilde '/vols_data/duncan_data'];

workingdir=[testdir]; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

clear spm_files_continuous spm_files_epoched;

% set up a list of SPM MEEG object file names (we only have one here)
spm_files_continuous{1}=[workingdir '/fMMspm8_EEG2_raw_0001.mat'];
spm_files_epoched{1}=[workingdir '/efMMspm8_EEG2_raw_0001.mat'];

% defining experimental conditions and contrasts to be calculated

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:


for ww=1:length(spm_files_epoched)

D=spm_eeg_load(spm_files_epoched{ww});

cnd=conditions(D);

mat=zeros(2,length(cnd));
for zz=1:length(cnd)
    str=cnd{zz};
    mat(1,zz)=str2num(str(5));   % here we assume that the format is LoadX_YY
    mat(2,zz)=str2num(str(7:end));
end

cat_var=unique(mat(1,:));

X=zeros(length(cnd),length(cat_var));

for k=1:length(cat_var)
    vect=find(mat(1,:)==cat_var(k));
    par_var=mat(2,vect);
    par_var=par_var-min(par_var)+1;
    X(vect,k)=par_var;
end

save([spm_files_epoched{ww} '.txt'],'X','-ascii');
design_matrix_summary={};
design_matrix_summary{1}=[spm_files_epoched{ww} '.txt'];



contrast={};

nc=length(cat_var);

for cont=1:nc
contrast{cont}=zeros(nc,1);
contrast{cont}(cont)=1; % main effects of a single condition (categorical variable)
contrast_name{cont}=['C' num2str(cont)];
end
for i=1:nc-1
    for j=i+1:nc
        cont=cont+1;
        contrast{cont}=zeros(nc,1);  % direct contrast between pairs of conditions
        contrast{cont}(i)=1; 
        contrast{cont}(j)=-1; 
        contrast_name{cont}=['C' num2str(i) '_vs_C' num2str(j)];
    end
end
for i=1:nc-1
    for j=i+1:nc
        cont=cont+1;
        contrast{cont}=zeros(nc,1);  % direct contrast between pairs of conditions
        contrast{cont}(i)=-1; 
        contrast{cont}(j)=1; 
        contrast_name{cont}=['C' num2str(j) '_vs_C' num2str(i)];
    end
end
cont=cont+1;
contrast{cont}=zeros(nc,1);  % average of all conditions
contrast{cont}(1:nc)=1; 
contrast_name{cont}='C_all';


%%%%%%%%%%%%%%%%%%%
%% DO SENSOR SPACE MULTI-BAND TIME-FREQ ANALYSIS USING OAT
% This fits the first-level GLM to sensor space data after it has been
% subject to a time-frequency transform

oat=[];
%oat.source_recon.D_continuous=spm_files_continuous{ww};
oat.source_recon.conditions=unique(cnd);
oat.source_recon.D_epoched{1}=spm_files_epoched{ww}; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[2 30]; % frequency range in Hz
oat.source_recon.time_range=[-0.1 0.5];
oat.source_recon.modalities={'EEG'};
oat.source_recon.method='none';
oat.source_recon.dirname=[oat.source_recon.D_epoched{1} '_hilbert_multiband'];
 

oat.first_level.tf_method='hilbert'; % can be morlet or hilbert
oat.first_level.tf_num_freqs=8; % we are keeping this unusally low in the practical for the sake of speed
oat.first_level.tf_hilbert_freq_res=4;
%oat.first_level.bc=[1 1 0];



oat.first_level.design_matrix_summary=design_matrix_summary;

% contrasts to be calculated:
oat.first_level.contrast={};
for i=1:length(contrast)
    oat.first_level.contrast{i}=contrast{i};
    oat.first_level.contrast_name{i}=contrast_name{i};
end

oat = osl_check_oat(oat);

%% now run the OAT

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

% load GLM result
stats_hilb=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% view the GLM design matrix (NOTE that column 1 is motorbikes, columns 2-4 are faces)
figure;imagesc(stats_hilb.x);title('GLM design matrix');xlabel('regressor no.');ylabel('trial no.');

%% visualise using Fieldtrip
% note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='EEG';
S2.first_level_contrast=[1];
S2.cfg.colorbar='yes';

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data_tf]=osl_stats_multiplotTFR(S2);

%% to do a topoplot averaged over 200 to 300 ms:
cfg.xlim        = [0.2 0.3];
cfg.ylim        = 'maxmin';
cfg.zlim        = 'maxmin';
cfg.interactive = 'no';

figure; ft_topoplotTFR(cfg,data_tf{1});




%%%%%%%%%%%%%%%%%%%
%% ANSWER/CHEAT: DO SENSOR SPACE SINGLE-BAND TIME-FREQ ANALYSIS USING OAT
% Here we zoom in on single frequency band (5 to 20 Hz) and look at the power of the
% acitivity over time in that band.
% To do this we change S2.num_freqs to be 1, and set the frequency range
% appropriately

oat=[];
% oat.source_recon.D_continuous=spm_files_continuous{ww};
oat.source_recon.conditions=unique(cnd);
oat.source_recon.D_epoched{1}=spm_files_epoched{ww}; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[5 20]; % frequency range in Hz
oat.source_recon.time_range=[-0.15 0.5];
oat.source_recon.modalities={'EEG'};
oat.source_recon.method='none';
oat.source_recon.dirname=[oat.source_recon.D_epoched{1} '_hilbert_singleband'];
 

oat.first_level.tf_method='hilbert'; % can be morlet or hilbert
oat.first_level.tf_num_freqs=1; % we are keeping this unusally low in the practical for the sake of speed
oat.first_level.tf_hilbert_freq_res=diff(oat.source_recon.freq_range);
oat.first_level.time_range=[-0.1 0.45];
oat.first_level.baseline_timespan=[-0.1 0];

oat.first_level.design_matrix_summary=design_matrix_summary;

% contrasts to be calculated:
oat.first_level.contrast={};
for i=1:length(contrast)
    oat.first_level.contrast{i}=contrast{i};
    oat.first_level.contrast_name{i}=contrast_name{i};
end



oat = osl_check_oat(oat);

%% now run the OAT

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

%% visualise using Fieldtrip
% note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='EEG'; 
S2.first_level_contrast=1;

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotER(S2);

end